Rare decays B — > X s (d)7 at the NLO 
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^ ■ Abstract 

O \ We present an independent calculation of 0(a s ) matrix elements of the QCD penguin operators 
& | for the decays B — ► X s r ( n'y. The Mellin-Barnes representation technique is applied for the 
^p" 1, systematic numerical evaluation of involved two-loop diagrams. The numerical effect of the 
calculated contributions is confirmed to be of order of 1%. We also present an updated numerical 
4^ analysis for B — > X s ^ d ^ decay rates. 
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1 Introduction 



The rare B decays play a special role in the phenomenology of weak decays. They are loop- 
induced in the standard model (SM) and thus are strongly affected by many contributions 
from the new sources of flavor violation in various extensions of the SM. The inclusive rare 
transitions are of particular interest as they are well approximated by the underlying partonic 
transitions. The corrections to this perturbative picture are calculable within the framework 
of Heavy Quark Effective Theory (HQET). 

Among the rare B-decays the inclusive B — > X s 7 transition is the most interesting one. 
Since the first observation of the decay by the CLEO [1] new measurements have been reported 
by ALEPH, BELLE and BABAR [2]- [5]. The CLEO still has the best determination [6]: 

B(B -> X sl ) = (3.21 ± 0.43 ± 0.27±g:io) x 10_4 - ( L1 ) 

while the current world average is [7] 

B(B -> X s7 ) = (3.40 ± 0.39) x 10" 4 . (1.2) 

The CKM-suppressed counterpart B — > X^ decay is also promising. However, no experimental 
observations of the b — > mediated transitions have been made so far. The upper experimental 
limits for the exclusive decays B — > p^ ^, B — > u^f [8] are higher than the corresponding 
theoretical predictions roughly by a factor of 2. 

As for the theoretical side, calculation of the NLO perturbative QCD corrections for the 
decay B — > X s ^ has been completed in [9]-[18]. These corrections lead to dramatic reduction of 
the renormalization scale dependence of the leading order (LO) results bringing the uncertainty 
in theoretical predictions for the decay rate to the level of 10% and opening the door for 
a rigorous comparison with experimental data. The leading order long distance corrections 
within the framework of HQET have been computed in [19]-[21]. Also the leading electroweak 
corrections are known [22]. The comparison of the current theoretical predictions with the 
available experimental data already presents strong constraints on the extensions of the SM 
(see e.g. [23] for a recent analysis). These analyses have been extended to the case of B — > X^ 
decay in [24]- [26]. 

Only about two years ago the last missing ingredients of the full NLO result- 0(a s ) ma- 
trix elements of the QCD-penguin operators, have been calculated in [18]. The calculation is 
complicated due to large number of contributing non-trivial two-loop diagrams, so a verifying 
re-calculation is desirable. In this paper we present an independent computation of these matrix 
elements. We apply the Mellin-Barnes representation technique for the systematic numerical 
calculation of the involved two-loop diagrams. This is in contrast to the method of [18] where 
the results are expressed in terms of a few integrals to be evaluated numerically. We confirm 
their numerical results and extend the analysis to the case of B — > X d ^ decay (see also [28] for 
a very recent similar analysis). 

The paper is organized as follows. In Section 2 we present the effective Hamiltonians for 
the decays B — > X s m^ and collect the final results for the matrix elements of the dimension-6 
operators. Section 3 is devoted to the details of our calculation method. In Section 4 we give 
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an update of numerical predictions for the branching fractions for B — > X s ^)1- Some of the 
useful formulae are collected in appendices. 

2 The decay rates for B — > X S ^)J 

The effective Hamiltonian for the decay b — > 57 (and the associated bremsstrahlung process 
6 — > reads 

W e //(6 - «7) = E CM - A« £ a (Of - OO) (2-1) 

where the CKM factors are given by 

Aj s) = V*V tb = -AX 2 (l - ^) (l - y + A 2 (P - , A« ee V^Kb = AA 4 (p-^).(2.2) 

The Wilson coefficients C(/x) are known to 0(a s ) precision. We refer to [16] and references 
therein for the explicit formulae. The dimension six effective operators can be chosen as 

O u x = (s Ll „T a u L ) (u L YT a b L ), O u 2 = (sl1»u l ) (u L ^b L ), 

Oi = {sLl»T a c L ){c L YT a b L ), 2 = (s Lllx c L ){c Ll »b L ), 

3 = (s L7 A) 4 = (s L 7„r°6 L ) E g (g7 M T a g), (2.3) 

5 = (s l7 m7,7p&l) 'EqiqYYYl), 6 = {s L l„lvl p T a b L ) E q (qrYYT a q), 

7 = Yi^m b (s L a^b R )F^, 8 = Y§^m b (s L a^b R )F^, 

The corresponding formulae for the decays b — > d'j(g) can be obtained from Eqs.(2.1),(2.3) with 
the obvious replacement s — > d The relevant CKM factors read 



A^ = V£V tt = AX 3 (l - (1-p + ifj), == V^Kft = AA 3 (p - 



(2.4) 



Taking into account that A is small (A ~ 0.22) one can safely neglect the terms proportional 
to in the case of b — > S7. In the case of b — > g?7 however these terms are important in 
particular being responsible for the generation of a large CP-asymmetry. We will keep those 
terms both in analytical and numerical formulae. 

The decay widths for the partonic transitions can be expressed as 

F[6 -> s ( d ^E 7 >E =(l-6)m b /2 = 



C 2 rv 
^ F^em 

32vr 4 
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s{d) 



' ^U^mN (\ DS(d) + e ewf + A'«>) , (2.5) 
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where we have separated the virtual and bremsstrahlung corrections. At the NLO we have 



D s( d ) = C (0)ef.f {flb) + 
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(2.7) 



The formulae for the bremsstrahlung coefficients /jj are presented in Appendix B. e ew represents 
the leading electroweak correction that can be extracted from [22]. 
The quantities r\ are given by 



n = - g r 2 , 

r 3 = 10.059 + |f in, 

r 5 = 185.841 + ^in, 



r 7 
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r2 = -ff> + 2(a(z)+6(z))-§*7r, 

r 4 = -1.017- ^j}i7r + 26(z), 

r 6 = -8.148 - ^§in + 12a(z) + 206(z), 

^8 



(2.8) 
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and r" = lim rj(i = 1,2). We will present the calculation of the matrix elements r$, ...,tq in 



some detail in Section 3. In Eq. (2.8) we have used the functions a(z),b(z) defined in [18] to 
encode z-dependent part of the matrix elements of the current-current operators 0\ and O2'. 
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where L = In z 1 . 

Following [23] we use the following formula for the branching ratios: 

The essential reason for introducing (2.11) is that the factor 



2 

c 



V ub Z T(B^ X c ev) 



Vr 



d> 



T(B -> X u ev) 



can be predicted with relatively small uncertainty using the T expansion [29]. We use C = 
0.575 x (1 ± 0.03) as estimated in [23]. As it follows from Eq. (2.11) 0(a s ) corrections for the 
inclusive decay B — > X u ev have to be included in our analysis. 

When doing the transition to B-meson decay rates in Eq. (2.11) from the corresponding 
partonic decay rates one has to include the leading 1/t?t| [19] and l/m 2 c corrections [21]. Note 
that in the case of B — > X^ additional long-distance contributions due to intermediate u- 
quark loops can have a sizeable effect. These contributions can not be reliably predicted but 
the model calculations indicate that they are not very large [30], [31]. In the present paper we 
will ignore them. 



3 Calculation of the virtual O (a s ) corrections for the op- 
erators O3, Oq. 

3.1 Organization of the calculation 

We start describing the virtual 0(a s ) contributions to the matrix elements of the operators 
3 ,...,Oq . The relevant Feynman diagrams are depicted in Fig.l. Note that the full set of 
contributing diagrams is larger, however, the calculation of the other diagrams can be avoided 
if one uses the gauge invariance of the result. Furthermore, the contributions of some of the 
diagrams are already included when one uses the effective Wilson coefficient Cj^ instead of 
C 7 . 

Note that for each of the diagrams there are two possible types of operator insertions as 
illustrated in Fig. 2. Apparently for type I insertions only b and s quark loops contribute while 
all five (u, d, s, c and b) quarks can give non-zero type II contributions. In our calculations we 
set the masses of the light quarks u, d, s to zero. 

Our calculation extends the method described in [17]. While for the light and c quark 
contributions it is straightforward, for the evaluations of the diagrams with b quark in the loop 
we use a different tactics to get numerical results. As we think that our method can be useful 
in other similar calculations we describe it in some detail in subsection 3.2. 

1 Originally these functions have been calculated in [17] with precision of 0(z 3 ) and later extended to include 
terms up to z 6 in [27]. We have verified these expressions using the method of [17] to include terms 0(z 4 ), 
0(z 5 ), and 0(z 6 ). 
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For a better organization of our computation we will use the following decomposition: 

3 = (s^Lb) 'Eiq-fLq) + {s llx Lb) Rq) = 3L + 3R , 

q 1 

4 = (s^LT a b) J2(qYLT a q) + (s 7 ^T a 6) ^(^KTq) = 4L + 4R 

q q 

and similarly for the operators O5, 0$. 

Following the ref. [17], as a first step in the calculation of the two-loop diagrams in Fig.l 
we evaluate the building blocks from Fig. 4. We present the full expressions of the one-loop 
building blocks in Appendix A. As we are going to illustrate the method of our calculation in 
the next subsection we will give here the expression for the contribution of 3 l with 6-quark 
loop to the building block Ip\ 

= -£,r( g )(l - e) ^V" / X)]1 ' € M - r^)L ± (3.1) 



where 

A = r 2 -m\jx(\ -x) +i5. (3.2) 



3.2 The method of the calculation 

We will illustrate our method of calculation on the example of diagram (1) from Fig.l for the 
type I insertion of the operator 3 l- Using (3.1) and applying the standard technique of the 
Feynman parameterization for the integration over the second 4-momentum we arrive to 2 

M(lc) = eQ f/^ F T(2e) exp 2me / dxdudvdy[x(l - x)] 1 ^' 1 ^ - v) e v x 

647T 4 J 

2 The details of the calculation of the same diagram with a charm-quark loop for the operator O2 are given 
in [17]. 
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Figure 2: Two types of operator insertions contributing to b — > s(d) r y. 



u(p') 



P 1 ^ + P 2 ^ + P 3 J— 

C 2e C 2e C l+2e 



u(p) 



(3.3) 



where the integration over the Feynman parameters x,y,u,v is over the interval [0;1] and C is 
given by 



C = m b v(l — v)u — m n —. + to 

q x[l — x) 



(3.4) 



The strategy in [17] (as well as in the similar calculations [32], [33]) is then to use the Mellin- 
Barnes representation for the integrals in eq. (3.3) which is given by 



(k 2 -M 2 Y (k 2 ) x T(\)2iri 



- ( ds(-M 2 /k 2 ) s T(-s)T(\ + s) 



(3.5) 



where A > and the integration path goes parallel to the imaginary axis in the complex s-plane 
hitting the real axis somewhere between —A and 0. In [17] the formula (3.5) is applied to the 
denominators in (3.3) using the identification 



k 2 <-> mlv(l — v)u ; M 2 <-> m 



(i -v)y 

x) 



2 

q x(l 



(3.6) 



Then the integrals over the Feynman parameters are expressed in terms of Euler's T-functions. 
To evaluate the integral over the Mellin-Barnes parameter s the integration path is closed in the 
right s-half-plane. When m 2 q /m 2 < 1/4 the integral over the half-circle vanishes so using the 
residue theorem the integral is expressed as the sum over residues in the integration contour. 
This in turn leads to a natural expansion in the (small) parameter m 2 q /m 2 . 

This procedure apparently won't work for the case q = b. The absence of a small expansion 
parameter seems to make the described technique unapplicable in this case. However we note 
that the appearance of the factor x(l—x) (0 < x < 1) in eq. (3.4) plays of a role of a suppression 
factor for the second term and it turns out that the 'opposite' identification 



M 



m 2 v{\ — v)u, k 2 <-> m 



x(l — x) 



(3.7) 
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does a better job in the case q = b. Just as in [17] then we first take the integrals over the 
Feynman parameters which after appropriate substitutions can be brought into the following 
forms: 



f 1 dww p , f 1 dww p (l-w) q = 
Jo Jo 



r(p+i)r( g + i) 

T(p + q + 2) 



Then to perform the integration over s we can close the integration path in the right half-plane 
and use the residue theorem. The relation x(l — x) < 1/4 guarantees then that the integral 
over the half-circle vanishes. 

So the integrals will be reduced to sums over the residues of T function that reside on the 
real axis (in complex s-plane): 

s = n, n — e, n — 2e (n — 0, 1, 2, ...), s = m/2 — 2e, (m = 3, 5, ...) (3.8) 

We note again that though there is no real expansion parameter, the terms with higher values 
of n (or m) are suppressed because of factor x{\ — x). So taking into account only first several 
values of n makes a good numerical approximation of the final result. 

It is easy to check the validity of our procedure for the diagram (1) from Fig. I. as it can be 
calculated analytically (see formula (5.6) of ref. [18])). Summing over the residues at n — 0, .., 6 
(see Eq.(3.8)) we will get real good approximation for the final result with the precision of 10~ 6 . 
We didn't observe such a fast convergence for the other diagrams, however the sum over the 
several first values of n in all of the cases gives a good approximation for the involved integrals. 
For each of the diagrams we will keep as many poles as it is necessary for getting a precision 
of at least 1(T 3 . 



3.3 Results of the calculation 

Using the method described above we are able to calculate the unrenormalized matrix element 
of the operator 3 . We find 

* = feK?r +9 - 387+ ^ <,r ) x s woH6> '~ <3 - 9) 

Looking at the expressions for the building blocks it is easy to see that the contributions from the 
operators can be expressed in terms of corresponding contributions from 3 with multiplicative 
factors that can be extracted from the formulae in Appendix A. The diagrams from Fig. 1.5 that 
do not contain the building blocks are easy to calculate without any numerical approximations. 
The results for the unrenormalized matrix elements of the remaining operators are 

M > = I -si; fir) " " 899 + 2 b{2) ~ wA x t^ ^- 



/5552 /m 6 \ 4e 1888. \ n s . ,,. 

Ms = UT7fif +68.688+— ^ x-MO^ 
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M 6 = (-H(y) 25 ' 956 + 12fl W + 206 W-W m ) x S (s7l ° 7|6)free 
where 

(S7|0 7 |6) tree = ^u(p')MRu(p). (3.10) 

These results are in a perfect agreement with those from [18]. 

The renormalization procedure is described in [18]. As we follow the same scheme we don't 
go into details. The final results for the renormalized matrix elements are given in Eq. (2.8). 



4 Numerical results 

In this section we will present the numerical results for the branching ratios of the decays 
B — > X s (d)7. In our numerical analysis we will use the following values taken from [35]: 

M w = 80.425 GeV, M z = 91.1876 GeV, m p t ole = 174.3 ± 5.1 GeV, a em = 1/137.036. (4.1) 

For the semileptonic branching ratio we will use the averaged value B(B — > X c £ + U£) exp = 
(10.90 ± 0.23)% from [36]. We use two-loop expression for the strong coupling constant a s (fi) 
with a s {M z ) = 0.1172 ±0.002 [35]. For the Wolfenstein parameters of the CKM matrix we will 
use [37] 

p = 0.162 ±0.046, r] = 0.347 ±0.027, A = 0.83 ± 0.02, A = 0.2240 ± 0.0036 (4.2) 

As it was pointed out in [23], in the phenomenological analysis of B — > X s ^ an important issue 
is the definition of the mass of charm quark m c in the matrix elements of the current-current 
operators Oi, 2 . The pole mass have been used in earlier works (e.g. [17]). However, as it was 
argued in [23], the use of the running mass m c (/i) with fj, varying between m c and rrib is more 
appropriate. We follow this choice using 

777 

— = 0.23 ±0.05 (4.3) 
m b 

in the matrix elements of the effective operators. Eq. (4.3) agrees with the numbers used in 
[28]. They are obtained using m c (m c ) = (1.25 ±0.10) GeV and m b = m\ s = (4.69 ±0.03) GeV 
[39]. The latter value for m b is used also anywhere else in the analysis. 

Using these input parameters and the photon energy cutoff E = 1.6GeV 3 we obtain for 
the branching ratios: 

B{B^X sl ) = (3.52l^(^)lH?W)^:^W H :^(CKM)±0.16(param)) xlO" 4 , (4.4) 
B{B^X dl ) = (l.34l^(^)^3W)i:^(^)l^ 8 7 (CKM)±0.07(param)) x 10" 5 . (4.5) 
3 This corresponds to d = 0.318 for m b = 4.69GeV. 
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Figure 3: The dependence of B(B — > X d ^y) on the CKM parameters p and rj. Thick 
(thin) lines correspond to the branching ratio with matrix elements of QCD pen- 
guin operators included (excluded). The dependence on p is shown for three values 
of rj: rj = 0.384 (upper curves), rj = 0.356 (middle curves) and rj = 0.328 (lower curves). 
Vertical lines illustrate the experimental bounds on the parameter p. 

The central values here correspond to the central values of the input parameters and /iw = 
Mw, /•*& = iTT-b- We allow pij, to vary between mb/2 and 2rrib and the matching scale piw between 
Mw and m\ ole . The last error sums up the uncertainties due to a s (M z ), m^ ole , B[B — > X c eu) exp 
and the factor C (see Eq. (2.11)). The effect of the matrix elements of the QCD penguin 
operators calculated in this paper is —1.0% for B — > X s ^ and —1.1% for B — > X d "f. These 
corrections lead to a slight increase of the renormalization scale dependence. Combining all the 
uncertainties we get 



As we have mentioned in Section 2 in our numerical calculations we have neglected the long- 
distance contributions due to presence of intermediate w-quark loops. Due to CKM-suppression 
these corrections are negligible for B — > X s ^y but can be sizeable for B — > X d ^. This is another 
source of uncertainty of theoretical prediction for B — > X d j that can be roughly estimated to 
be of order of 10% (see e.g. [24] for a more detailed discussion). 

The results of our numerical analysis are somewhat lower than those from [23], [18]. Note 
that the central values in Eq. (4.5) will decrease by about 1% if we use the values from [23] for 
the input parameters. The numerical discrepancy is explained by the difference in methods used 



B(B - X sl ) 
B(B - X dl ) 



(3.45 ±0.30) x 10 
(1.27 ±0.24) x 10 



(4.6) 
(4.7) 
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in the analysis. The main difference is that in [23] the charm-quark and top-quark contributions 
are split and each of them is treated separately. 

As we can see from Eq. (4.6) the uncertainty due to CKM parameters while being negligible 
for B — > X s 7 is still dominating for B — > X<f). In Fig. 3 we illustrate the dependence of 
B(B — > X(fi) on p and rj. 

To conclude, we have presented an independent calculation of the 0(a s ) corrections to the 
matrix elements of the QCD penguin operators O3, ...,Oq. The effect of these corrections is 
about 1% decrease for the branching ratios. The updated predictions for the branching ratios 
are (3.45 ± 0.30) x 10~ 4 for B -> X s7 and (1.27 ± 0.24) x 10~ 5 for B -> X d7 . 
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Appendix A: Formulae for one-loop building blocks 

In this appendix we will give the complete formulae of the one-loop building blocks for all 
different contributions described in section 3.1. All the calculations are done using the dimen- 
sional regularization in d = 4 — 2e dimensions using the NDR scheme with fully anticommuting 
75. For both diagrams index (3 is to be contracted to the gluon propagator and a to the photon 
polarization vector, r and q denote the 4-momenta of the gluon and photon respectively. 




Figure 4: One-loop building blocks used in our calculation. 

For type I insertion the contribution of the g-quark loop for the operators On and Om 
(q — s,b; i = 3, 6) are respectively 

I$(L) = -ClDl{L)^Y{e){l - e)^e^ J " ^ (rrf - r^)L ±, (A.l) 



I$(R) = -C\D\{R)^Y{e)e^e^) xdx[x ^ x)] ~ ' m q (f^ - lp f)R{q) * (A.2) 



where R{s) = L, R{b) = R, 

A = r 2 - m 2 q /x(l -x)+ iS. (A.3) 
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The operator-dependent color coefficients Cj and Dirac algebra induced factors A are given by 

Cl = C{ = \, Ci = C^ = -l/(2N c ); 
Di(L,R)=D{(L,R) = l, 

Dl(D = Dq(L) = 4(4 - e - e 2 ), Di(R) = Dq(R) = 4(5 - 3c - e 2 ). 
For the type II insertions Ip(R) apparently vanishes so we have (q = u, d, s, c,b; i = 4, 6) 



rq(H) _ MH)(T\ 



± 

/ 



dx[x(l — x)] 1 



(rpf 



(A.4) 



where 

Cf = C{ = 0, C? = C» = (1/2); 
D 4 7/ = 2, £> 6 " = 4(5-3e). 

Similarly, for the type I contributions of operators OjL and for the building block J a p 
we have 



16tt 2 



£(a, ft r)A* 5 + £(a, (3, q)Ai & - E((3, r, q)—Ai 23 

qr 



-E(a, r, q)^Ai 25 - E(a, r, q)^Ai 26 



qr 



qr 



A 



(A.5) 



= ^ m<?e IQn^ ^ is dxdyC l € ( 8 (9ap(qr) - r a qp) 



x 



(dxyBYiR) - 2Bf{R)) + 2(2 + e)Bf (R)(fwJ + i lal pf) 



where 



E(a, P,r) = - {icclpf - flpla) ■ 



h\ (A.6) 



(A.7) 



The operator-dependent color coefficients Ai and Dirac algebra induced factors Bi are given by 

Ai = Al = l, A{ = Ai = -l/(2N c )- 
Bi I (L,R)=Bi I (L,R) = l, j = 1,2,3, 

Bi(L) = Bq(L) = 4(4 - 56 - e 2 ), Bl\R) = B?(R) = 4(1 + c - e 2 ), 
Bf(R) = Bf(R) = 4(1 - 3e - 3e 2 ), Bf(R) = Bf{R) = 4(1 - 7e - e 2 ). 

Finally the quantities Aj are defined as in [17] 

Ai 23 = -Ai 26 = 8(gr) f cfody [xye(l + e)r(e) exp^e)// 6 *?- 1 " 

J s 

Ai 25 = -8(qr) f dxdy - x)e(l + e)r(e) exp(7 £ e)// 2e C~ 1_e 



(A.8) 



(A.9) 



Ai 5 = Ai 23 , Ai 6 



(qr) 



Ai 25 + Ai 26 , 
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where C is given by 

C = m 2 q — 2xy(qr) — x(l — x)r 2 — i5 
In the case of type II contributions we have 



M H hr\ — AH nII (T \ e 9sQq 



E{a, /3, r)At 5 + E(a, (3, q)Ai & - E((3, r, q)—Ai 23 

qr 



-E{a, r, q) 1 - l At 25 - E(a, r, q)^Ai 26 
qr qr 



L-, 

2' 



where 



(A.10) 



Bi 1 
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= Bi l = Bi l =^ Bq — —24/ (6 — d) 
= 1/2- 



Appendix B: Formulae for bremsstrahlung corrections 

Here we list the formulae for the bremsstrahlung coefficients fff 4 5 : 
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Voz 

Az 2 



5 



(l-S)/z 



dt (1 - zt) 



mi 

t 2 



9 



(l-5)/z 



t 



dt ( G(t) + - ) 



+ 



l/z 



1/z 
(1-8)/. 



dt (1 - zt)' 



G(t) + 1 



t 



(l-5)/z 



dt(l-zt) (G(t) + - 



4 
9 
1 
3 
1 

27 



25 3 

105 + 5 2 - — + 5(5 -4) In 5 



21n^(5 2 + 25 + 41n(l-5)) 



+4Li 2 (l-5)- 
5 



2tt 



- 5(2 + 5) In 5 + 81n(l - 5) - -5 3 + 35 2 + 75 
3 3 



8z 
27 



dt(l - zt) ' 



r 



(l-5)/ 2 



dt(l - zf) 



2 , G*(t) 



r 



-i-5 (3 - 35 + 5 2 ) 



i22 — 



27 
2 

"81 



5(5 2 -3), 



fB.ll 



4 Note that our definitions of quantities ft? slightly differ from the standard definitions of (see e.g. [16]) 
5 The formula for the coefficient /gg is the only place where the mass of the final state quark has to be retained 
due to logarithmic divergence. This colinear divergences can be resummed to all orders in perturbation theory 
[38]. However, for realistic values of m s and md the contribution of fss is under 1%, so we will ignore this issue 
using mt,/m s = 50 and nib /ma = 1000 in our analysis. 
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where the function G(t) is defined in Eq. (40) of [16]. The other coefficients can be expressed 
in terms of those defined above: 



for the quantities defined in Eqs. (B.l), (B.2). In the bremsstrahlung corrections we have used 
the same value m c /m b = 0.23 ± 0.05 as in Eq. (4.3) . 
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